 
                     
            disp('simulating')
            
            disp('Reminder: now simulating for coaldev,  exact')
            disp([coaldev exact])
            disp('now simulating for parameters')
            disp([Qbeta,Qrho])
            disp('vss')
            disp(vss)
            
            
            Npanel=Nvill*(vss);
            csim0=nan(Nvill,vss,T);
            income0=nan(Nvill,vss,T);
            income0ind=nan(Nvill,vss,T);
            phisim0=nan(Nvill,vss,T);
            Yvill=nan(1,T);
            Yagg=nan(Nvill,T);
            Ssim=nan(Nvill,vss,T);
            gammarfun=nan(Nvill,vss,T);
            phinew=nan(Nvill,vss,T);
            wwnew=nan(Nvill,vss,T);
            indSjj=nan(T,vss);
            gammaraux=nan(Nvill,vss,T);
            t=1;
             
            % initial values of income and consumption (set =
            % income)
            rng(vss,'twister');
            shockinit=randi(length(incomevals),Nvill,vss);
            income0ind(:,:,1)=shockinit;
            
            income0(:,:,1)=income(income0ind(:,:,1));
            
            csim0(:,:,1)=income0(:,:,1);
            income=incomevals;
            
            Yagg(:,t)=sum(income0(:,:,t),2);
            
        tic   
            IDIO_DIST_sim=[];             
while t<T
                t=t+1;
                %%%%simulate next period income taking into account
                %%%% persistence across time
                rng(t,'twister');
                shock=rand(Nvill,vss);
                for k=1:agent
                    if length(incomevals)==3
                    rr=(find(shock(:,k)<=cumQ(income0ind(:,k,t-1),1)));
                    income0ind(rr,k,t)=1;
                    rr=(find(cumQ(income0ind(:,k,t-1),1)<shock(:,k) & cumQ(income0ind(:,k,t-1),2)>=shock(:,k)));
                    income0ind(rr,k,t)=2;
                    rr=(find(cumQ(income0ind(:,k,t-1),2)<shock(:,k) & cumQ(income0ind(:,k,t-1),3)>=shock(:,k)));
                    income0ind(rr,k,t)=3;
                    income0(:,k,t)=income(income0ind(:,k,t));
                    elseif length(incomevals)==2
                    rr=(find(shock(:,k)<=cumQ(income0ind(:,k,t-1),1)));
                    income0ind(rr,k,t)=1;
                    rr=(find(cumQ(income0ind(:,k,t-1),1)<shock(:,k) & cumQ(income0ind(:,k,t-1),2)>=shock(:,k)));
                    income0ind(rr,k,t)=2;
                    income0(:,k,t)=income(income0ind(:,k,t));    
                    end
                end
              
                Yagg(:,t)=sum(income0(:,:,t),2);
                gammarfun(:,:,t)=csim0(:,:,t-1)./(csim0(:,1,t-1)*ones(1,agent));
                
                
                
                for ivill=1:Nvill
                %%%find the state we're in
                if vss==2
                check=[SS(:,1)==income0ind(ivill,1,t)].*[SS(:,2)==income0ind(ivill,2,t)];
                elseif vss==3
                check=[SS(:,1)==income0ind(ivill,1,t)].*[SS(:,2)==income0ind(ivill,2,t)].*[SS(:,3)==income0ind(ivill,3,t)];
                elseif vss==4
                check=[SS(:,1)==income0ind(ivill,1,t)].*[SS(:,2)==income0ind(ivill,2,t)].*[SS(:,3)==income0ind(ivill,3,t)].*[SS(:,4)==income0ind(ivill,4,t)];
                end
                jsim(ivill,t)=find(check==1);  
                if vss==2
                    for g=1:2
                phinew(ivill,g,t)=interp1(gammagrid(:,4), phisol(:,jsim(ivill,t),g), gammarfun(ivill,2:agent,t));
                   end
                elseif vss==3
                phinew(ivill,:,t)=interpolationphi(phi,jsim(ivill,t), gammarfun(ivill,2:agent,t), vss);
                
                elseif vss==4
                phinew(ivill,:,t)=interpolationphi_4agents(phi,jsim(ivill,t), gammarfun(ivill,2:agent,t), vss);
                end
                
                end
                
                
                
                          
                gammaraux(:,:,t)=((1+phinew(:,:,t))./(1+phinew(:,1,t)*ones(1,vss))).^(1/rho).*csim0(:,:,t-1)./(csim0(:,1,t-1)*ones(1,vss));
                if vss==3
                for ivill=1:Nvill
                wwnew(ivill,:,t)=interpolationphi(wsb,jsim(ivill,t), gammaraux(ivill,2:agent,t), vss);
                end
                end
                
                
                csim0(:,:,t)=(gammaraux(:,:,t)./(sum(gammaraux(:,:,t),2)*ones(1,vss))).*(Yagg(:,t)*ones(1,vss));
                
                
end
            csim0keep=csim0;
            income0keep=income0;
            AA=[];
            BB=[];
            CC=[];
            DD=[];
            for t=1:T
                AA=[AA reshape(csim0(:,:,t)',Nvill*(vss),1)];
                BB=[BB reshape(income0(:,:,t)',Nvill*(vss),1)];
                CC=[CC reshape(phinew(:,:,t)',Nvill*(vss),1)];
                DD=[DD reshape(wwnew(:,:,t)',Nvill*(vss),1)];
            end
            csim0=AA;
            income0=BB;
            phinew0=CC;
            wwnew0=DD;
            